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Abstract 

Screening and surveillance are routinely used in medicine for early detection of disease and 
close monitoring of progression. Biomarkers are one of the primarily tools used for these 
tasks, but their successful translation to clinical practice is closely linked to their ability 
to accurately predict clinical endpoints during follow-up. Motivated by a study of patients 
who received a human tissue valve in the aortic position, in this work we are interested in 
optimizing and personalizing screening intervals for longitudinal biomarker measurements. 
Our aim in this paper is twofold: First, to appropriately select the model to use at time 
t, the time point the patient was still event-free, and second, based on this model to select 
the optimal time point u > t to plan the next measurement. To achieve these two goals 
we develop measures based on information theory quantities that assess the information we 
gain for the conditional survival process given the history of the subject that includes both 
baseline information and his/her accumulated longitudinal measurements. 
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1 Introduction 


Decision making in medicine has become increasingly complex for patients and practitioners. 
This has resulted from factors such as the shift away from physician authority toward shared 
decision making, unhltered information on the Internet, new technology providing additional 
data, numerous treatment options with associated risks and benehts, and results from new 
clinical studies. Within this context medical screening procedures are routinely performed 
for several diseases. Most prominent examples can be found in cancer research, where, 
for example, mammographic screening is performed for the detection of breast cancer in 
asymptomatic patients, and prostate-specihc antigen (PSA) levels are used to monitor the 
progression of prostate cancer in men who had already been diagnosed with the disease. In 
general, the aim of screening procedures is to optimize the benehts, i.e., early detection of 
disease or deterioration of the condition of a patient, while also balancing the respective 
costs. In this paper we are interested in optimizing screening intervals for asymptomatic or 
symptomatic patients that are followed-up prospectively (the latter possibly after a medical 
procedure). The setting we focus on is better explained via Figure 1. More specihcally, this 
hgure depicts a hypothetical patient who has been followed-up up to time point t, and has 
provided a set of longitudinal biomarker measurements. Using this information, physicians 
are interested in events occurring within a medically relevant time interval {t,t + Af]. If 
the probability of surviving beyond t + At is lower than the pre-specihed constant c{t), 
then the physician is supposed to immediately take action to improve the prospects for this 
patient. However, when the chance of surviving until t + At is higher than c{t), then the 
physician may not initiate treatment or require an invasive test, and instead wishes to gain 
a better understanding of the disease progression using additional biomarker measurements. 
The decision to request an additional measurement requires careful weighing of the potential 
benehts and costs to the patient and health care system. Hence, we are interested in carefully 
planning the timing of the next measurement in order to maximize the information we can 
gain on disease progression, and at the same time minimize costs and patient burden. 

Statistical methods for optimizing medical screening strategies have hourished in the past 
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Figure 1: Illustration of a patient’s longitudinal profile up to time t, and conditional survival 
function for whom we wish to plan an extra measurement. Physicians are interested in 
events occurring within the time interval {t, t + At]. Constant c{t) defines the threshold 
above which the physician requires extra information before making a decision. 


years. These methods have been primarily based on Markov and multistate models (Son- 
nenberg and Beck 1993; Parmigiani 1993; Lee and Zelen 1998; Parmigiani 1998, 1999, 2002; 
Schaefer et ah 2004; Alagoz et ah 2010). The appealing feature of these models is that they 
can simultaneously handle different medical options and select the one that provides the 
most gains at an acceptable cost (cost-effectiveness). Here we will instead work under the 
framework of joint models for longitudinal and time-to-event data (Henderson et ah 2000; 
Tsiatis and Davidian 2004; Rizopoulos 2012). An advantageous feature of joint models is 
that they utilize random effects and therefore have an inherent subject-specific nature. This 
allows one to better tailor decisions to individual patients, i.e., personalize screening, rather 
than using the same screening rule for all. In addition, contrary to the majority of screen- 
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ing approaches that are based on a dichotomized version of the last available biomarker 
measurement to indicate the need for a further medical procedure (e.g., PSA greater than 
4.0 ng/mL), joint models utilize the whole longitudinal history to reach a decision without 
requiring dichotomization of the longitudinal outcome. Our methodological developments 
couple joint models with information theory quantities, following similar arguments as in 
Commenges et ah (2012). In particular, our aim in this paper is twofold: First, to appropri¬ 
ately select the joint model to use at time f, i.e. the time point the subject of interest is still 
event-free, and second, based on this model to select the optimal time point u > t to plan 
the next measurement. We measure optimality by the amount of information gained for 
the survival process given the history of the subject of interest that includes both baseline 
information and his/her accumulated longitudinal measurements. Depending on the charac¬ 
teristics of the disease under study, the physician may request that u lies within the interval 
{t,t -|- At] or set another upper time limit after which he/she is not willing to wait before 
obtaining the next measurement. 

The motivation for this research comes from an aortic valve study conducted by the 
Department of Cardio-Thoracic Surgery of the Erasmus Medical Center in the Netherlands. 
This study includes 285 patients who received a human tissue valve in the aortic position in 
the hospital from 1987 until 2008 (Bekkers et al. 2011). Aortic allograft implantation has 
been widely used for a variety of aortic valve or aortic root diseases. Major advantages as¬ 
cribed to allografts are the excellent hemodynamic characteristics as a valve substitute; the 
low rate of thrombo-embolic complications, and, therefore, absence of the need for anticoag¬ 
ulant treatment and the resistance to endocarditis. A major disadvantage of using human 
tissue valves, however is the susceptibility to degeneration and the concomitant need for 
re-interventions. The durability of a cryopreserved aortic allograft is age-dependent, leading 
to a high lifetime risk of re-operation, especially for young patients. Re-operations on the 
aortic root are complex, with substantial re-operative risks, and mortality rates in the range 
of 4-12%. A timely well-planned re-operation has lower mortality risks compared to urgent 
and re-operative surgery. It is therefore of great interest for cardiologists and cardio-thoracic 
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surgeons to plan in the best possible manner echocardiographic assessments that will inform 
them about the future prospect of a patient with a human tissue valve in order to optimize 
medical care, carefully plan re-operation and minimize valve-related morbidity and mortality. 

The rest of the paper is organized as follows. Section 2 briefly presents the joint modeling 
framework that sets the basis for our later methodological developments. Section 3 presents 
the information theory quantities that can be used to select the joint model at a particular 
follow-up time, and Section 4 the corresponding quantities for planning the next longitudinal 
measurement. Section 5 illustrates the use of these information measures in the Aortic Valve 
study, and Section 6 concludes the paper. 

2 Joint Model Specification 

In this section we present a general dehnition of the framework of joint models for longitu¬ 
dinal and survival data that will be used later on for planning the optimal visit schedule. 
Let Vn = = l,...,n} denote a sample from the target population, where T* 

denotes the true event time for the Ath subject, Ci the censoring time, Ti = min(T*, Ci) the 
corresponding observed event time, and = I{T* < Ci) the event indicator, with /(■) being 
the indicator function that takes the value 1 when T* < Ci, and 0 otherwise. In addition, 
we let Ui denote the n* x 1 longitudinal response vector for the f-th subject, with element 
yu denoting the value of the longitudinal outcome taken at time point tu, I = 1,..., n^. 

To accommodate different types of longitudinal responses in a unihed framework, we 
postulate a generalized linear mixed effects model. In particular, the conditional distribution 
of yi given a vector of random effects hj is assumed to be a member of the exponential family, 
with linear predictor given by 

k [E{yi{t) I bi}] = r]i{t) = xj{t)f3 + zj{t)bi, (1) 

where k{-) denotes a known one-to-one monotonic link function, and yi{t) denotes the value of 
the longitudinal outcome for the f-th subject at time point t, Xi(t) and Zi(t) denote the time- 
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dependent design vectors for the fixed-effects f3 and for the random effects hj, respectively. 
The random effects are assnmed to follow a mnltivariate normal distribntion with mean zero 
and variance-covariance matrix D. For the snrvival process, we assnme that the risk for an 
event depends on a fnnction of the snbject-specihc linear predictor r]i{t) and/or the random 
effects. More specihcally, we have 


hi{t \T-ii{t)^Wi) = lim Pr{t < T* < f-|-At I T* > t,'Hj(t), 

At—>0 

= ho{t) exp[y'^Wi + f{r]i{t),bi,cx}], t > 0, (2) 

where 'Hi(t) = {? 7 j(s ),0 < s < t} denotes the history of the nnderlying longitndinal process 
np to t, ho(-) denotes the baseline hazard fnnction, Wi is a vector of baseline covariates with 
corresponding regression coefficients 7 . Fnnction /(■), parameterized by vector ck, specihes 
which components/featnres of the longitndinal ontcome process are inclnded in the linear 
predictor of the relative risk model. Some examples, motivated by the literatnre (Brown 
2009; Rizoponlos and Ghosh 2011; Rizoponlos 2012; Taylor et ah 2013; Rizoponlos et ah 
2014), are: 




f{'Hi{t),bi,cx} 


f{'Hi{t),bi,cx} 


f{'Hi{t),bi,cx} 


ar]i{t), 


air]i{t) + a 2 r][{t), with 77 '(t) 

a / T]i{s)ds, 

Jo 

oJbi. 


drjiit) 

dt 


These formnlations of /(■) postnlate that the hazard of an event at time t may be associated 
with the nnderlying level of the biomarker at the same time point, the slope of the longitn¬ 
dinal prohle at t, the accnmnlated longitudinal process up to t, or the random effects alone. 
Finally, the baseline hazard function ho(-) is modeled flexibly using a B-splines approach. 


i.e.. 


Q 

\ogho{t) = -fhofi + ^-fho,qBq{t,v), 

9=1 


( 3 ) 
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where Bq{t,v) denotes the g-th basis function of a B-spline with knots vi,... ,vq and 
the vector of spline coefficients. To avoid the task of choosing the appropriate number 
and position of the knots, we include a relatively high number of knots (e.g., 15 to 20) 
and appropriately penalize the B-spline regression coefficients for smoothness using the 
differences penalty (Eilers and Marx 1996). 

For the estimation of joint model’s parameters we use a Bayesian approach based on 
Markov chain Monte Carlo (MCMC) algorithms. The expression for the posterior distribu¬ 
tion of the model parameters given the observed data is derived under the assumptions that 
given the random effects, both the longitudinal and event time process are assumed inde¬ 
pendent, and the longitudinal responses of each subject are assumed independent. Formally 
we have. 


p{yi,Ti,5i 

\b^,e) 

= Piyi\bi,0)p{Ti,6i\bi,e), 

( 4 ) 

PiVi 

\b^,e) 


( 5 ) 


i 


where 0 denotes the full parameter vector, and p(-) denotes an appropriate probability 
density function. Under these assumptions the posterior distribution is analogous to: 

n Tii 

p{e,b)^ nn piUii I bi, 0) p{Ti, 5i I hi, 0) p{bi I 6) p{6), (6) 

i=l 1=1 


where 


piPii I bi,0) = exp 


ViiAiibi) 


c{ipii{bi)} /a{p)-d{yii, 



with '^ii(hi) and p denoting the natural and dispersion parameters in the exponential family, 
respectively, c(-), a(-), and d{-) are known functions specifying the member of the exponential 
family, and for the survival part 


p(Ti, 6i I bi, 6) 


hi{Ti I 'Hi(Ti, bi)Y^ exp 


r-Ti 


hi(s I 'hLi 


s, b. 


)) ds 
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with hi{ ) given by (2). The integral in the definition of the survival function 


Si{t I bi 


Wi 


= exp< — 


ho(s)exp + f{r]i{s),bi, 0 '}]ds 


( 7 ) 


does not have a closed-form solution, and thus a numerical method must be employed for its 
evaluation. Standard options are the Gauss-Kronrod and Gauss-Legendre quadrature rule. 

The penalized version of the B-spline approximation to the baseline hazard can be fitted 
by specifying for the improper prior (Lang and Brezger 2004): 


Phho I 'Th) 


oc 



where Th is the smoothing parameter that takes a Gamma(l, r/i^) prior distribution, with 
a hyper-prior Ths ~ Gamma(10“^, 10“^), which ensures a proper posterior distribution for 
(Jullion and Lambert 2007), K = AjAj. + 10“®I, with denoting the r-th difference 
penalty matrix, and p{K) denotes the rank of K. 


3 Choosing the Model to Use at Time t 

As motivated in Section 1, our aim is to decide the appropriate time point to plan the next 
measurement for a new subject j from the same population (i.e., with j ^ {i = 1,... ,n}), 
whose survival chance at a particular time point t lies within the interval [c{t), 1] (see Fig¬ 
ure 1). Recent research in the field of joint models (Yu et ah 2008; Proust-Lima and Taylor 
2009; Rizopoulos 2011, 2012; Taylor et al. 2013) has shown how such conditional survival 
probabilities can be computed in order to identify for which subjects we need to proceed 
in planning the next measurement. In particular, the fact that biomarker measurements 
have been recorded up to t, implies that subject j was event-free up to this time point, and 
therefore we are interested in: 

7lj{u I t) = Pt{T* > u\T* > 


(8) 


where yj{t) = {yjitji)] 0 < tji < t,l = 1,..., rij} denotes the longitudinal measurements of 
subject j taken before t, and Vn = {Tj, 6i, y^;i = 1,... ,n} denotes the original dataset. 

As we have seen in Section 2, in general, defining a joint model entails choosing an 
appropriate model for the longitudinal outcome (i.e., baseline covariates and functional form 
of the time effect), an appropriate model for the time-to-event (i.e., baseline covariates), and 
how to link the two outcomes using function /(■). Hence, a relevant question is which model 
to use for deciding when to plan the next measurement of subject j. Standard approaches for 
model selection within the Bayesian framework include the deviance information criterion 
(DIG) and (pseudo-) Bayes factors. However, a potential problem with these methods is that 
they provide an overall assessment of a model’s predictive ability, whereas we are interested 
in the model that best predicts future events given the fact that subject j was event-free up to 
time point t. To identify the model to use at t we follow similar arguments as in Commenges 
et ah (2012), and focus on the conditional density function of the survival outcome given the 
longitudinal responses and survival up to t. More specihcally, we let Ai = {Mi,..., Mk} 
denote a set of K joint models htted to the original dataset Vn- Due to the fact that the 
choice of the optimal model in Ai is based on a hnite sample, an important issue we need 
to address is overfitting. Hence, to obtain an objective estimate of the predictive ability of 
each model we will use the cross-validatory posterior predictive conditional density of the 
survival outcome, which for model M^ is defined as p{T* \ T* > t,yi(t),Vn\i, , where 
Dn\i = {Til, 5ii,yii;i' = 1,... ,i — l,i + 1,... ,n} denotes the version of the dataset that 
excludes the responses of the Ath subject. Similarly we can define the same distribution 
under model M*, which denotes the model under which the data have been generated (M* 
is not required to be in Ad). Using conventional quantities of information theory (Cover and 
Thomas 1991), we select the model Mk in the set Ad that minimizes the cross-entropy: 


CEk{t) = U -log p{T* I T* > t,yyt),V^\i,Mk} 


where the expectation is taken with respect to p(T* \ T* > t,yi(t),Vn\i, M*). An estimate 
of CEfc(f) that accounts for censoring can be obtained using the available information in the 
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sample at hand. We term this estimate the cross-validated Dynamic Conditional Likelihood'. 


1 X ^ 

cvDCLfc(t) = — > t)\ogp{Ti,5i \ Tt > t, W(t), (9) 

1 

2=1 

where nt = > ^)- turn, an estimate of p{Ti,6i \ T* > t,yi{t),'Dn\i, Mk] can be 

obtained from the MCMC sample of model Mk by utilizing a similar computation as for the 
conditional predictive ordinate statistic (Chen et ah 2008, condition on Mk is dropped in 
the following expressions but is assumed): 


p{Ti,6i I Ti > t,yi{t),Vn\i] 


p{Vn\i,Ti > t, W(t)} 
p{Ti, 6i,Ti> t, yi{t),Vn\i} 


f P{T^n\i 

6)p{Ti > t^yyt) 


1" p{Ti > t,yi{t) 

1 

\e)p{e) 

1 


de j 


0} p{Vn) 


p{Ti,6i I Ti > t,yi{t),6} 


p{6 I Vn) de, 


( 10 ) 


where 6 includes here both the parameters of the joint model 6, as dehned in Section 2, and 
the random effects. Hence, based on the MCMC sample from [6 \ D^] and identity (10) we 
can obtain a Monte Carlo estimate of (9) as 


- 1 

cvDCL(f) = -Vj(T, >f) log 


G 


G 

T 


9 = 




( 11 ) 


where 6^^'^ denotes here the g-th. realization from the MCMC sample = 1,... ,G}. 

Calculation of p{Tj, 6i \ Ti > f, W(t), is based on the conditional independence assump¬ 
tions (4) and (5). Due to space limitations and the fact that similar calculations will be 
required in the following section we do not show them in detail here. 
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4 Scheduling Patterns based on Information Theory 


4.1 Planning of the Next Measurement 

Having chosen the appropriate model to use at time t, our aim next is to decide when to 
plan the following measurement of subject j. Similarly to the previous section, the quantity 
that we will use for making this decision is the posterior predictive conditional distribution 
p{T* I T* > t,yj{t),Vn}. This distribution again has a cross-validatory flavor because 
j ^ {i = I,... ,n}. We let yj{u) denote the measurement of the longitudinal outcome at 
the future time point u > t. For the selection of the optimal u we would like to maximize 
the information we gain by measuring yj{u) at this time point, provided that the patient 
was still event-free up to u. If the event occurs before time u for subject j, then there 
is no gain in information. Using concepts from Bayesian optimal designs (Verdinelli and 
Kadane 1992; Clyde and Chaloner 1996), this can be formally expressed using the utility 
function (conditioning on baseline covariates Wj is assumed in the following expressions but 
is dropped for notational simplicity): 


U(m I f) = £'< Ai log 


Pj'^j I {yjit),yjiu)},'Dn) 

p{T* I T* > u,yj{t),Vn} 


+ A2/(t; > u) 


( 12 ) 


where the expectation is taken with respect to the joint predictive distribution [T*,yj{u) \ 
T* > t,yj(t),Vn\. Making the connection to the information theory concepts we used 
in Section 3, the hrst term in (12) is the expected Kullback-Leibler divergence between 
the posterior predictive conditional distributions with and without this extra measurement, 
namely 


EKL(m \t)=E 


Y 


E. 


'T*\Y 


log 


pj'^j I {yj{t),yj{u)},Vn) 
p{T* I T* > u,yj{t),Vn} 


(13) 


log 


p{T* I T; > u, {yj{t),yj{u)},Vn) ^ I ^ m ^ ^ t . 7 ^* 

p{T- 1 t ; > u,yM,v„} I > i. {yAt),m(n)},T>„) dT, 


y<p{yj{u) 


T* > t,yj{t),Vn}dyj{u) 
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Heuristically, for u close to t, in general, we would expect minimal information gain. Anal¬ 
ogously, as u moves further away from t we would expect to gain information, but the value 
of u that maximizes this information gain will depend on the characteristics of the subject’s 
projected longitudinal profile. However, the ‘cost’ of waiting up to time u, is the fact that 
we also inevitably increase the risk that the patient will experience the event. Because we 
would not like to wait up to a point that it would be late for the physician to intervene (e.g., 
the event already happened), we need to counterbalance these two aspects, which motivates 
including the second term in (12). The expectation with respect to this second term is the 
conditional survival probability 7ij{u \ t) introduced in (8). Hence, (12) becomes 

U(m I t) = AiEKL(m I t) + X2'iTj{u I t). 

The purpose of the nonnegative constants Ai and A 2 is to weigh the contribution of the 
cumulative risk as opposed to the information gain, and to take into account that these two 
terms have different units. In practical applications, elicitation of these constants for trading 
information units with probabilities can be difficult. However, using the equivalence between 
compound and constrained optimal designs (Cook and Wong 1994; Clyde and Chaloner 
1996), it can be shown that for any Ai and A 2 , there exists a constant k G [0,1] for which 
maximization of (12) is equivalent to maximization of EKL(m | t) subject to the constraint 
that 'n'j{u \ t) > K. Under this equivalent formulation elicitation of k is relatively easier; 
for example, a logical choice is to set k equal to the constant c(t), introduced in Figure 1, 
for deciding whether the physician is supposed to intervene or not. In turn, a choice of 
c{t) can be based on medical grounds or on time-dependent receiver operating characteristic 
(ROC) analysis. Then, the optimal value of u can be found by maximizing EKL(m | t) in 
the interval (t,t“^], where = {u : 7ij{u \ t) = k}. In practice, it also may be reasonable 
to assume that there is an upper limit of the time interval the physician is willing to 
wait to obtain the next measurement, in which case = minju : 7ij{u \t) = K, The 

advantageous feature of this approach is that the optimal u is found by conditioning on both 
the baseline characteristics Wj and accumulated longitudinal information 3^j(t) of subject j. 
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This allows one to better tailor screening to the subject’s individual characteristics compared 
to approaches that only utilize the baseline covariates and the last available longitudinal 
measurement (Markov assumption). 

4.2 Estimation 

Estimation of (12) proceeds by suitably utilizing the conditional independence assump¬ 
tions (4) and (5). More specihcally, the density of the joint distribution [T*,yj{u) \ T* > 
t,yjit),Vn] can be written as the product of the predictive distributions, 

p{t; I T;>t,{y,{t),y,{u)},Vn) = j p{t; \ t; >t,{y,{t),y,{u)}, 0 )pio \ v^)de, 


and, 

pivjiu) I T* > t,yj{t),Vn} = j p{yj{u) I T* > t,yj{t),e}p{6 \ Vn)do. 

Analogously, the hrst term in the right-hand side of each (4.2) and (4.2) can be written as 
expectations over the corresponding posterior distributions of the random effects of subject 
J, i-e., 


p{yj{u) I T* > t,yj{t),e} 


jp{yj{u) I hj,e}p{hj 


T* > t,yj(t),0} dbj, 


and. 


p{T;\T;>t,{y,{t),yj{u)},e) 

= jP{Tj I T* > t, bj, 6)p{bj I T* > t, {yj{t),yj{u)},6} dbj, 
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respectively. In a similar manner, the density of the predictive distribution [T* \ T* > 
u, {yj{t),yj{u)],Vn] is written as: 


p{T;\T;>u,{y,{t),y,{u)}M 

= j j I > u,bj,e)p{bj I T* > u, {yj{t),yj{u)},6)p{6 \ Vn)dbjdO, 


with 


p(t; 


T’ > u, bj, e) 


1 'Hj{t,bj,0))ex.p^ 

1-. 

j hj (s Hj (t, bj , 0)) ds 1 

exp j 

[-J 

pu 

1 h,{ 

0 

s 

Hj(t, bj, 0)) dsj 



and the hazard function is given by (2), in which we have explicitly noted that the history 
of the longitudinal process 'Hj(t,bj,6^ is a function of both the random effects and the 
parameters, because it is estimated from the mixed model (1). By combining the above 
equations we can construct a Monte Carlo simulation scheme to obtain an estimate of (12) 
for a specihc value of u. This scheme contains the following steps that should be repeated 
q = 1,... ,Q times: 

~ ^ (i) 

Step 1: Simulate 0, 0, and {6 = independently from the posterior distribution 

of the parameters [6 \ V„]. 

Step 2: Simulate bj from the posterior of the random effects [bj \ T* > t,yj{t),e]. 

Step 3: Simulate yj{u) from [yj{u) \ bj,6]. 

Step 4: Simulate bj and {bj ,i = 1,... ,L} independently from the posterior of the random 

^ (i) 

effects [bj \ T* > t,{yj{t),yj{u)},6] and [bj \ T* > u,{yj{t),yj{u)},6 ], respectively. 

Step 5: Simulate T* from [T* \ T* > u, bj, 0]. 
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L 

Step 6: If f* > u, compute | t) = log|L“^ , where 

1=1 

4,'> = hj{f; I nAf;,bfA))expl^-£’ h,{s \ n,(s,bf,eA) *}. 

and 

A^[^ = exp I “ y I ('5’ ^ I ’ 

otherwise set EKL*''^^(m | t) = 0. 

The sample mean over the Monte Carlo samples provides the estimate of EKL(m | t), i.e., 

EKL(m I t) = EKL^'^^(m I t). The posterior distribntion of the random effects in 

~ ■■ -p) 

Steps 2 and 4 is not of a known form, and therefore the realizations hj, bj and {bj ,i = 
1,..., L} are obtained nsing a Metropolis-Hastings algorithm, with proposal distribntion a 
mnltivariate Stndent’s-t distribntion with mean b = arg max^{logp(h | T* > t,yj{t),d)}, 
variance-covariance matrix V = [—d^ \ogp{b \ T* > t,yj(t),0)/db~^db\^^^^ and fonr 
degrees of freedom, where 6 denotes the posterior means based on [6 \ Vn]- Step 3 just 
entails simulating a realization from the mixed effects model. In Step 5 again the distribution 
[T* I T* > t,bj,6] is not of standard form, and hence we simulate the realization T* using 
the inversion method, i.e., hrst we simulate n ~ f/(0,1), and then using a line-search method 
we hnd the T* for which Sj{T*, 'bj, '6)/Sj{t, 'bj, '6) = v, with Sj{-) given by (7). Finally, in 
Step 6 it is sufficient to compute the numerator of (13) to obtain the optimal u. An estimate 
of 7^j{u I t) can be obtained from a separate Monte Carlo scheme (see e.g., Rizopoulos 2011) 

hj[s I ds 

hj[s I 'Hj[s,b^f\6^^^)) ds 

with {0^^\g = 1,...,G} a sample from [6 \ P„], and b^f^ a sample from [bj \ T* > 
t,y,it),e^% Based on the Monte Carlo sample of EKL(m | t) a 95% credible interval 
van be constructed using the corresponding Monte Carlo sample percentiles. Using this sim- 
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ulation scheme we take as optimal u the value u = argmax^g^^^^ | t)}, subject 

to the constraint that 9j{u \t)>K, with ui,, Umax denoting a grid of values in (t, t^^]. 

5 Analysis of the Aortic Valve Dataset 

We return to the Aortic Valve dataset introduced in Section 1. Our aim is to use the existing 
data to build joint models that can be utilized for planning the visiting patterns of future 
patients from the same population. In our study, a total of 77 (27%) patients received a 
sub-coronary implantation (SI) and the remaining 208 patients a root replacement (RR). 
These patients were followed prospectively over time with annual telephone interviews and 
biennial standardized echocardiographic assessment of valve function until July 8, 2010. 
Echo examinations were scheduled at 6 months and 1 year postoperatively, and biennially 
thereafter, and at each examination, echocardiographic measurements were taken. By the 
end of follow-up, 1262 assessments have been recorded, with an average of 4.3 measurements 
per patient (s.d. 2.4 measurements), 59 (20.7%) patients had died, and 73 (25.6%) patients 
required a re-operation on the allograft. Here we are interested in the composite event 
re-operation or death, which was observed for 125 (43.9%) patients. 

We start by dehning two sets of joint models that we will use for planning measurements 
of future patients. Each set considers a separate longitudinal outcome that measures the 
functioning of the aortic valve, namely, the square root transform of aortic gradient (con¬ 
tinuous outcome) and the aortic regurgitation (binary outcome). Preliminary analysis using 
plots of the subject-specihc trajectories of the two outcomes, and htting linear mixed effects 
and mixed effects logistic regression models suggested nonlinear prohles for aortic gradient 
and linear evolutions for the logit (logit (x) = log{a;/(l — x)}) of the probabilities of aortic 
regurgitation. Hence, we let ^ denote the continuous aortic gradient longitudinal outcome. 
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and 7/j 2 the binary aortic regurgitation outcome and postulate the models 




(/^0,l + ^iO,l) + 1) + (/^2,1 + ^j2,l)-Bn(t, 2) 

+ /?3,iAge. + /34,iFemalei + ei(t), 


where Bn{t, {1, 2}) denotes the B-spline basis for a natural cubic spline with boundary knots 
at baseline and 19 years and one internal knot placed at 3.7 years (i.e., the median of the 
observed follow-up times), Female denotes the dummy variable for females, and 

logit [Pr{|/j^ 2 (t) = 1 }] = Wo ,2 + bio, 2 ) + Wi ,2 + biiWt + / 32 , 2 Age- + / 33 ^ 2 Feni 3 -len 


Both hj 1 and hj ^2 are assumed to follow multivariate normal distributions with mean zero 
and covariance matrices Di and D 2 , respectively. For the survival process we consider four 
relative risk models, each positing a different association structure between the two processes, 
namely: 


Ml,k ■ 

hi,kit) 

to 

hi,k it) 

M3,k ■ 

hi,kit) 

MA,k '■ 

hi,kit) 

M5,k ■ 

hi,kit) 


hoWt) exp{7i,fcAge. -f 72,fcFemalei -F ai^kVi,k{t)}, 
hoWt) exp{7i,fcAge. -F 72 ,A:Femalei -F a2,kVi,k(t)}y 
ho,kW) exp{7i,fcAge. -F 72,fcFemalei -F ai^kr]i,kW) + a2,kr]'i,kW)]y 


ho,kit) exp<^ 7 i,fcAge. -F 72 ,fcFemalei -F ai,k 



ho,kit) exp(7i,feAge. -F 72,fcFemalej -F cxjbi,k ), 


where k = 1,2 denotes the relative risk models corresponding to the two sets of joint models 
for aortic gradient and aortic regurgitation, respectively. The baseline hazard was approxi¬ 
mated using penalized B-splines approach described in Section 2. For all parameters we took 
standard prior distributions (Ibrahim et ah 2001; Brown and Ibrahim 2003; Brown 2009). 
In particular, for the vector of hxed effects of the longitudinal submodel (3, the regression 
parameters of the survival model 7 , the vector of spline coefficients for the baseline hazard 
and for the association parameter a we used independent univariate diffuse normal 
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priors. For the variance of the error terms cx^ we take an inverse-Gamma prior, while for 
covariance matrices we assnmed an inverse Wishart prior. 

Tables 1 and 2 present the posterior means and 95% credible intervals of the parameters of 
the joint models for aortic gradient and aortic regurgitation fitted to the Aortic Valve data. 
With respect to the parameters of the longitudinal component we do not observe great 



Mi,l 

Value (95% Cl) 

M2,1 

Value (95% Cl) 

^3,1 

Value (95% Cl) 

M4,i 

Value (95% Cl) 

M5,1 

Value (95% Cl) 

Intercept 

B-splnl 

B-spln2 

Age 

Female 

(j 

dll 

d2l 

dsi 

d,22 

d,32 

d.33 

3.67 (3.33; 4.00) 
3.28 (2.79; 3.83) 
2.70 (2.33; 3.13) 
-0.02 (-0.02; -0.01) 
0.17 (-0.06; 0.40) 
0.61 (0.58; 0.65) 
0.64 (0.48; 0.83) 
-0.69 (-1.44; -0.06) 
-0.43 (-1.05; 0.09) 
13.61 (9.66; 18.74) 
8.14 (4.40; 13.62) 
7.04 (3.11; 12.99) 

3.66 (3.32; 4.01) 
3.37 (2.84; 3.89) 
2.77 (2.34; 3.21) 
-0.02 (-0.02; -0.01) 
0.17 (-0.06; 0.39) 
0.62 (0.58; 0.65) 
0.62 (0.47; 0.80) 
-0.51 (-1.19; 0.10) 
-0.40 (-0.97; 0.15) 
13.87 (9.79; 19.14) 
9.46 (5.84; 14.20) 
8.08 (4.62; 13.25) 

3.68 (3.34; 4.03) 
3.31 (2.81; 3.84) 
2.72 (2.35; 3.13) 
-0.02 (-0.02; -0.01) 
0.17 (-0.04; 0.41) 
0.62 (0.58; 0.65) 
0.63 (0.48; 0.84) 
-0.64 (-1.38; -0.04) 
-0.43 (-1.00; 0.13) 
13.50 (9.76; 18.14) 
8.05 (4.55; 12.82) 
6.63 (3.32; 11.71) 

3.67 (3.34; 4.02) 
3.26 (2.76; 3.75) 

2.68 (2.28; 3.08) 
-0.02 (-0.02; -0.01) 

0.18 (-0.05; 0.41) 
0.61 (0.58; 0.65) 
0.64 (0.47; 0.82) 
-0.69 (-1.42; -0.05) 
-0.46 (-1.10; 0.11) 
13.91 (9.85; 19.09) 
8.70 (4.61; 13.66) 
7.80 (3.55; 14.06) 

3.67 (3.32; 4.02) 
3.34 (2.85; 3.83) 
2.72 (2.31; 3.15) 
-0.02 (-0.02; -0.01) 
0.18 (-0.04; 0.41) 
0.62 (0.58; 0.66) 
0.61 (0.46; 0.81) 
-0.53 (-1.23; 0.12) 
-0.34 (-0.96; 0.28) 
13.06 (9.26; 17.83) 
8.01 (4.09; 13.47) 
7.31 (3.06; 15.57) 


Mi_2 

^2,2 

^3,2 

^4,2 

^5,2 


Value (95% Cl) 

Value (95% Cl) 

Value (95% Cl) 

Value (95% Cl) 

Value (95% Cl) 

Intercept 

-2.01 (-3.74; -0.34) 

-2.03 (-3.80; -0.41) 

-1.89 (-3.43; -0.32) 

-2.03 (-3.64; -0.43) 

-1.98 (-3.57; -0.32) 

time 

0.38 (0.28; 0.48) 

0.39 (0.29; 0.49) 

0.39 (0.29; 0.50) 

0.37 (0.27; 0.47) 

0.40 (0.29; 0.49) 

Age 

-0.03 (-0.07; -0.00) 

-0.03 (-0.07; 0.00) 

-0.03 (-0.07; -0.00) 

-0.03 (-0.07; -0.00) 

-0.03 (-0.07; -0.00) 

Female 

1.40 (0.29; 2.55) 

1.40 (0.28; 2.56) 

1.26 (0.23; 2.24) 

1.43 (0.31; 2.51) 

1.35 (0.25; 2.44) 

dll 

14.29 (8.77; 22.97) 

14.04 (9.16; 21.57) 

11.01 (7.19; 16.10) 

13.39 (8.18; 20.15) 

12.64 (8.08; 18.83) 

d,2i 

-0.56 (-1.47; 0.17) 

-0.38 (-1.11; 0.24) 

0.06 (-0.58; 0.63) 

-0.48 (-1.31; 0.22) 

-0.18 (-0.92; 0.46) 

d22 

0.47 (0.27; 0.80) 

0.44 (0.24; 0.73) 

0.48 (0.27; 0.84) 

0.52 (0.30; 0.90) 

0.43 (0.25; 0.74) 


Table 1: Estimated coefficients and 95% credibility intervals for the parameters of the lon¬ 
gitudinal submodels fitted to the Aortic Valve dataset. The top part of the table refers to 
the results for aortic gradient and the bottom part for aortic regurgitation, dij denotes the 
ij-th element of the corresponding covariance matrix of the random effects. 



Mi.i 

Value (95% Cl) 

M 2.1 

Value (95% Cl) 

M3,i 

Value (95% Cl) 

M4,i 

Value (95% Cl) 

Ms.i 

Value (95% Cl) 

Age 

Female 

Ql 

a2 

0^3 

0.02 (0.01; 0.04) 
-0.09 (-0.48; 0.29) 
0.19 (0.09; 0.30) 

0.02 (0.01; 0.04) 
-0.02 (-0.43; 0.36) 

2.40 (1.07; 3.82) 

0.02 (0.01; 0.04) 
-0.09 (-0.49; 0.30) 
0.16 (0.02; 0.30) 
0.78 (-1.76; 2.99) 

0.02 (0.00; 0.03) 
-0.10 (-0.51; 0.29) 
0.02 (0.01; 0.03) 

0.02 (0.01; 0.04) 
0.01 (-0.43; 0.45) 
-0.03 (-0.59; 0.51) 
0.43 (0.26; 0.65) 
-0.02 (-0.75; 0.81) 


Value (95% Cl) 

^2.2 

Value (95% Cl) 

^ 3,2 

Value (95% Cl) 

^4,2 

Value (95% Cl) 

^ 5,2 

Value (95% Cl) 

Age 

Female 

Ql 

0:2 

0.01 (-0.00; 0.03) 
-0.02 (-0.43; 0.36) 
0.03 %.00; 0.07) 

0.01 (-0.00; 0.03) 
-0.01 (-0.44; 0.40) 

1.13 (0.32; 2.32) 

0.01 (-0.00; 0.03) 
0.07 (-0.40; 0.52) 
-0.09 (-0.18; -0.01) 
2.31 (0.84; 4.17) 

0.01 (-0.00; 0.02) 

0.00 P0.40; 0.41) 
0.001 (-0.002; 0.004) 

0.02 (0.00; 0.03) 
-0.02 (-0.51; 0.43) 
-0.06 %.22; 0.04) 
1.46 (0.43; 3.63) 


Table 2: Estimated coefficients and 95% credibility intervals for the parameters of the sur¬ 
vival submodels htted to the Aortic Valve dataset. The top part of the table refers to the 
results for aortic gradient and the bottom part for aortic regurgitation. 


differences between the htted joint models for both aortic gradient and aortic regurgitation. 
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For the relative risk submodels and for aortic gradient we see that when the terms rii^i{t) 
and 77 ' i{t) are independently included in the model both of them show to be related with the 
risk of re-operation/death, however, when both of them are included in the model, the slope 
term 77 ' does not seem to be associated with the hazard of the event anymore. Moreover, 
the area under the subject-specific longitudinal profile up to time t is also strongly related 
to the risk for event at same time point, indicating that not only features of the trajectory 
at time t are important but also at previous time points. For aortic regurgitation we observe 
that when combined, both the current value and the current slope term are associated with 
the risk of the composite event, whereas past values (i.e., the cumulative effect) do not seem 
to contribute in predicting future events. 

We continue by investigating the performance of the fitted joint models with respect to 
predicting future events after different follow-up times. Table 3 shows the DIG and the value 
of QYDCL(f) X Ut for f = 5, 7, 9,11 and 13. From the set of joint models fitted with aortic 




Mi.i 

M2,1 

A73,i 

M44 

M5.1 

DIG 


7237.26 

7186.18 

7195.57 

7268.45 

7186.34 

QvDCL(t = 

5) X nt 

-387.67 

-348.96 

-380.82 

-441.69 

-356.31 

cvDCL(t = 

7) X nt 

-342.75 

-309.69 

-336.95 

-391.78 

-315.26 

cvDCL(t = 

9) X nt 

-289.11 

-260.44 

-284.13 

-332.06 

-264.95 

c;vDCL(t = 

11) X nt 

-233.56 

-208.11 

-229.27 

-270.28 

-212.06 

(.vDCL(t = 

13) X nt 

-177.10 

-156.88 

-173.58 

-206.40 

-160.03 



Mi 2 

-^2,2 

A73.2 

^4,2 

^ 5,2 

DIG 


4166.76 

4114.27 

4049.99 

4186.72 

4061.53 

(;vDGL(t = 

5 ) X nt 

-370.46 

-362.97 

-328.10 

-358.95 

-365.31 

(;vDGL(t = 

7) X nt 

-328.84 

-321.25 

-288.47 

-318.74 

-322.70 

CTDGL(t = 

9) X nt 

-278.23 

-272.55 

-244.73 

-269.81 

-273.23 

CTDGL(t = 

11) X nt 

-223.45 

-218.67 

-195.18 

-216.21 

-219.19 

cvDGL(t = 

13) X nt 

-170.24 

-167.14 

-147.95 

-164.21 

-167.95 


Table 3: DIG and QyDGL(t) evaluated at hve time points for the fitted joint models. 


gradient as the biomarker we observe that models M 2 ,i and have almost identical DIG 
values but according to QvDGL(t) model M 24 provides the best predictions for future events. 
For aortic regurgitation model M 3_2 is selected as the best model according to both DIG and 
QyDGL(t). However, as a side note we observe that model has a lower DIG value than 
M 4^2 but the (;;yDGL(t) was higher for model M 4^2 for all t, which again suggests that the 
DIG, as a global measure, may not always select the model that more accurately predicts 
future events after a particular time point t. 
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To illustrate the use of the methodology developed in Section 4 for selecting the next 
screening time point, we focus on the continuous marker and two patients from the Aortic 
Valve dataset, namely Patient 81 who is a male of 39.8 years old, and Patient 7 who is a 
male of 46.8 years old and neither had an event. These patients have been explicitly selected 
because of the shape of their longitudinal trajectories of the square root aortic gradient, 
shown in Figure 2. More specihcally. Patient 7 exhibits a slightly increasing trajectory of 



Time (years) 


Figure 2: Longitudinal trajectories of the square root aortic gradient for Patients 7 and 81. 


aortic gradient indicating that he remained in relatively good condition during the follow-up, 
whereas Patient 81 showed an increasing prohle of aortic gradient indicative of a deterioration 
of his condition. We should note that in accordance with a realistic use of (12) for planning 
the next measurement, the joint models presented above have been htted in the version of 
the Aortic Valve dataset that excluded those two patients. Tables 4 and 5 show the values 
of EKL(m I t) (i.e., the hrst term of (12)), and of 7ij{u \ t) (i.e., the second term of (12)) for 
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these two patients computed in a dynamic manner, namely, after each of their longitudinal 
measurements (i.e., time point t). The upper limit of the interval was set as the 

minimum of = 5 years and the time point u for which 7ij{u \ t) was equal to k = 0.8. 
Then the optimal u was taken as the point that maximizes EKL(m | f) in a grid of hve 
equidistant values in the interval (f, All calculations have been based on joint model 
M 2 ^i. For Patient 7, who remained relatively stable during the hrst period of his follow-up. 


t 

-t 

u 

EKL{u 1 t) 

TTj(u 1 t) 


0.3 

5 

1.3 

-2.488 

0.989 




2.3 

-2.205 

0.977 




3.3 

-2.227 

0.963 




4.3 

-2.197 

0.948 




5.3 

-2.153 

0.933 

* 

1.3 

5 

2.3 

-2.176 

0.987 




3.3 

-2.127 

0.973 




4.3 

-1.963 

0.957 




5.3 

-2.088 

0.940 




6.3 

-1.868 

0.922 

* 

3.3 

5 

4.3 

-2.143 

0.978 




5.3 

-1.436 

0.956 

* 



6.3 

-2.015 

0.932 




7.3 

-1.913 

0.905 




8.3 

-1.634 

0.875 


5.3 

3.9 

6.1 

-2.143 

0.968 




6.9 

-1.603 

0.935 

* 



7.6 

-1.879 

0.900 




8.4 

-1.967 

0.864 




9.2 

-1.750 

0.827 


7.1 

2.1 

7.5 

-1.888 

0.969 

* 



7.9 

-1.993 

0.937 




8.4 

-1.963 

0.906 




00 

bo 

-1.901 

0.875 




9.2 

-2.104 

0.844 


10.6 

1.7 

11.0 

-0.243 

0.970 




11.3 

-0.156 

0.940 

* 



11.6 

-1.397 

0.910 




12.0 

-0.542 

0.880 




12.3 

-2.184 

0.849 



Table 4: The values of EKL(m | t) (i.e., the hrst term of (12)), 7ij{u \ t) (i.e., the second 
term of (12)), for Patient 81 computed in a dynamic manner, after each of his longitudinal 
measurements (i.e., time point f), and for u set as a regular sequence from f to f-|-3 of length 
hve. The last column denotes which time point is selected to plan the next measurement. 


we observe that EKL(m | t) always suggests waiting the full hve years up to measurement 
seven; however, after the seventh measurement and because his prohle started to slightly 
increase and his relative older age, the length of the interval (f, for hnding the optimal 
u decreases, which in turn suggests to wait after each visit less time before requesting the 
next echocardiography. For Patient 81 on the other hand we observe that for the hrst two 
visits, where he exhibited relatively low levels of aortic gradient, EKL(m | t) suggests to wait 
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for a longer period years before requesting the next measurement. In accordance, is set 
to five years because the probability that he will survive more than five years was higher 
than K = 0.8. After the third measurement however, when it starts to become clear that the 
condition of this patient deteriorates, EKL(m | t) recommends to wait less time to request 
the next echocardiography. Accordingly from the fourth measurement onwards the length 
of the interval decreases and in the last measurement is set to 1.7 years that also 

indicates the deterioration of the patient. 

6 Discussion 

In this paper we have presented two measures based on information theory that can be used 
to dynamically select models in time and optimally schedule longitudinal biomarker mea¬ 
surements. Contrary to standard screening procedures that require all patients to adhere to 
the same screening intervals, the combination of the framework of joint models for longitu¬ 
dinal and time-to-event data with these two measures allows one to tailor screening to the 
needs of individual patients and dynamically adapt during follow-up. 

Motivated by the aim of the Aortic Valve dataset, in this work we have concentrated 
on a single endpoint. Nevertheless, using recent advances in the joint modeling literature 
(Andrinopoulou et ah 2014; Dantan et al. 2011) the same ideas can be relatively easily 
extended to multiple endpoints (competing risks) or intermediate states combined with one 
or several endpoints (multistate settings). In such settings the first term of (12) will remain 
the same, but the second term 7ij{u \ t) will need to be replaced by the corresponding 
cumulative incidence function. Furthermore, our approach also could be combined with 
concepts, such as the value of information (Koerkamp et al. 2006) to incorporate the relative 
costs of measuring the biomarker. 

Finally, to facilitate the computation of QvDCL(f) and EKL has been implemented in 
functions cvDCLO and dynInfoO, respectively, available in package the JMbayes (ver¬ 
sion 0.7-0; Rizopoulos 2015) for the R programming language (freely available from the 
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Comprehensive R Archive Network at http://cran.r-project. org/package=JMbayes). 

A working example of how these functions should be used can be found in the supplemen¬ 
tary material. 
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t 

t^P -1 

u 

EKL(u 1 t) 

TTj{u 1 t) 


3.5 

5 

4.5 

- 2.363 

0.983 




5.5 

- 2.354 

0.964 




6.5 

- 2.306 

0.943 




7.5 

- 2.279 

0.919 




8.5 

- 2.065 

0.892 

* 

4 

5 

5.0 

- 2.366 

0.982 




6.0 

- 2.422 

0.963 




7.0 

- 2.210 

0.941 




8.0 

- 2.214 

0.915 




9.0 

- 2.106 

0.886 

* 

5 

5 

6.0 

- 2.416 

0.981 




7.0 

- 2.470 

0.959 




8.0 

- 2.261 

0.933 




9.0 

- 2.300 

0.904 




10.0 

- 2.146 

0.872 

* 

5.7 

5 

6.7 

- 2.455 

0.979 




7.7 

- 2.410 

0.956 




8.7 

- 2.311 

0.928 




9.7 

- 2.352 

0.897 




10.7 

- 2.207 

0.861 

* 

6.7 

5 

7.7 

- 2.426 

0.976 




8.7 

- 2.389 

0.949 




9.7 

- 2.328 

0.917 




10.7 

- 2.304 

0.882 




11.7 

- 2.235 

0.841 

* 

7.7 

5 

8.7 

- 2.403 

0.973 




9.7 

- 2.394 

0.941 




10.7 

- 2.307 

0.906 




11.7 

- 2.272 

0.865 




12.7 

- 2.171 

0.819 

* 

9.1 

4.3 

9.9 

- 2.355 

0.971 




10.8 

- 2.335 

0.939 




11.7 

- 2.201 

0.903 




12.5 

- 2.130 

0.863 




13.4 

- 2.096 

0.817 

* 

11.3 

3.1 

11.9 

- 2.202 

0.969 




12.5 

- 2.044 

0.935 




13.1 

- 2.003 

0.897 




13.8 

- 1.908 

0.856 

* 



14.4 

- 1.920 

0.810 


13.5 

1.9 

13.9 

- 1.919 

0.968 




14.3 

- 1.855 

0.933 




14.7 

- 1.742 

0.897 




15.1 

- 1.655 

0.858 

* 



15.4 

- 1.719 

0.817 


15.1 

1.4 

15.4 

- 1.769 

0.968 




15.7 

- 1.731 

0.934 




15.9 

- 1.632 

0.899 




16.2 

- 1.560 

0.862 




16.5 

- 1.536 

0.824 

* 


Table 5: The values of EKL(m | t) (i.e., the hrst term of (12)), Tij{u \ t) (i.e., the second 
term of (12)), for Patient 7 computed in a dynamic manner, after each of his longitudinal 
measurements (i.e., time point t), and for u set as a regular sequence from f to f+ 3 of length 
five. The last column denotes which time point is selected to plan the next measurement. 
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